Modeling of zonal electrophoresis in plane channel of complex shape 
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Abstract 

The zonal electrophoresis in the channels of complex forms is considered mathematically with 
the use of computations. We show that for plane S-type rectangular channels stagnation regions 
can appear that cause the strong variations of the spatial distribution of an admixture. Besides, 
the shape of an admixture zone is strongly influenced by the effects of electromigration and by 
a convective mixing. Taking into account the zone spreading caused by electromigration, the 
influence of vertex points of cannel walls, and convection would explain the results of electrophoretic 
experiments, which are difficult to understand otherwise. 
Keywords: Microchip, electrophoresis. 

Introduction 

During the last ten years or so the intense use of various microchips aimed to separate mix- 
tures by an externally imposed electric field, to control micromixing and chemical reactions 
has been flourished (see e.g. [|l], |, |, |, |, i 0, g |, |ll 0, [l^, |li |3 



The industrial use of microchannels for the efficient separation of mixtures is well known as 
the technology called Lab-on-a-Chip. The most effective control of mass transfer in micro- 
fabricated fluid devices can be achieved with the use of electrokinetic phenomena such as 
electrophoresis and electroosmosis. The crucially important part of related research is com- 
puter modeling that helps to improve the design of microchips, to understand the processes 
involved, and to enhance experimental methods. 

The modeling of the electrophoretic separation of a mixture represents a challenging 
problem due to the large number of physical phenomena involved into the mass transfer 



driven by an electric field. One can count here such phenomena as diffusion, chemical re- 
actions, dependence of electrical conductivity on concentrations, electroosmosis. Joule heat, 
convection, etc. It should be noticed that many papers devoted to the transport phenom- 
ena in microchannels take into account only diffusion, electroosmosis, and the Taylor-Aris 
dispersion. As the result these papers leave out of account some essential nonlinear effects 
that appear due to the dependence of electrical conductivity on component concentrations; 
it is well known that these effects significantly change zone shapes |2l|, H, |3l |2|, |2|, || 



and even can trigger a so called substance-lock effect |^ (for an experimental verification 
see p8[| ). 

The effects of electric field singularities that occur near the vertex points of electrophoretic 
chamber walls still have not been understood and mathematically described. For a simple 
case of a plane cross-shaped channel this singularity is of the order 0{r~^^^), where r is a dis- 
tance from the vertex point of the reflex angle 37r/2. The related distortions of zone shapes 
are described in [29 and experimentally justified in |^|. In addition, transport processes can 
be effected by convective mixing, which can drastically deform the final stage of a separation 



process. The role of convection in electrophoresis is described in [|31| , |32| , |33| , p4| , pSj p6| . It 



is apparent that one can weaken convection by the choosing of an appropriate orientation 
of a microchip in the gravity field. Nevertheless, taking convection into account can not be 
avoided for high precision experiments. 

This paper is devoted to the computer modeling that reveals the effects of such key 
factors as vertex points of channel walls, electromigration, and convection on zone distortions. 
We present only the results of a small part of our numerical experiments that have been 
carried out with the use of a specially created interactive program for the modeling of zonal 
electrophoresis. The main result is the identifying of the parameter intervals when the 
distortions of moving zones are the most significant. These data can be very useful for the 
planning of new experiments and for the designing of electrophoretic chambers. 



1 Mathematical Model 

The mathematical models of electrophoresis (and in particular zonal electrophoresis) are 
well known |^, ^ |3^. The dimensionless governing equations describing both the 
motion of separated (by the action of an electric field) substances and fluid convection (in 
the Oberbeck-Boussinesq approximation) are: 

^ = -Vp + /iAv - fc^/3fcCfc, divv = 0, (1) 

k=l 

^ + divzfc = 0, ik = -£|7fc|Vcfc + -^kCkE, (2) 



E = -Vip, j = crE, div(crVv9) = 0, (divj = 0), 



(3) 




> 0, 



(4) 



Here v and p are velocity and pressure, E and ip are the strength and the potential of an 
electric field, is the k-th concentration [k = 1, . . . ,n), — the fiux of concentration, 
j — the density of electric current, a — mixture conductivity, (Tq — the conductivity in 
the absence of admixtures (the conductivity of a buffer solution), /i — fluid viscosity, 7^ — 
electrophoretic mobility, l3k — the coefficient that appears in the linear dependence of density 
on concentration, — the coefficient that appears in the linear dependence of conductivity 
on concentration, e — characteristic diffusion coefficient {e\'-fk\ are diffusion coefficients), 
k — the unit vector of the z-axis that is anti-parallel to the gravity. The dimensionless 
variables used are described in 



In the presented model the effects of Joule heat and electroosmosis have been neglected. 
The conductivity of a mixture has been modeled by the expression that is natural for 
zonal electrophoresis; it is accepted that a mixture contains components with constant con- 
centrations that represent so called buffer solution (such that at = its conductivity is 
(To). The concentrations of separated substances (samples) are assumed to be small enough. 
We should emphasize that the coefficients can have different signs (positive or negative); 
from a physical viewpoint ak < means that this particular substance has lower specific 
conductivity than that of a buffer solution. When this substance enters into a solution, it 
'replaces' the buffer substances and the conductivity of a mixture is decreasing (for details 
see 



21, 24, 25 



Let us consider the domain shown in Fig. |1 . 1| . We accept that its boundary is rigid with 
non-leak conditions for a liquid and for the concentrations 



V 



0, ik ■ n = 0, k = 1, 



n. 



(5) 



where n is the unit vector to the boundary. 
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Figure 1.1: Electrophoretic chamber 



The potential (f is prescribed on the parts BC and GF, while the rest of the boundary 
is insulated: 



BC 



GF 



(6) 



On 



0. 



(7) 



CDEF,BAHG 

At the initial instant a fluid is still and the initial distributions of admixtures are prescribed: 

c^(a:, z), k = 1 n. 



t=0 



0, Cfc 



(8) 



2 Qualitative Analysis of the Problem 

There are at least four factors that influence the distortion of a zone as well as admixture 
concentrations. The first one is diffusion that causes the spreading of electrophoretic zones. 
For small concentrations (the case of analytical electrophoresis) nonlinear effects are weak 
and diffusion makes the separation of admixtures difficult; its influence is well studied and 
described in almost all handbooks on zonal electrophoresis (see [21, 0). The second factor 
is the electromigration spreading of zones. It reveals itself for high concentrations (for 
example for preparative electrophoresis) and described in details in |j2l], ^ Recall, 
that in one- dimensional case the evolution of an initially rectangular concentration profile 
for a single admixture follows the patterns shown in Fig. where two upper pictures give 
initial distributions for positive and negative ai, while two bottom pictures present some 
later stages of their developments. 
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Figure 2.1: The zone evolution 



The analysis of diffusionless {e = 0) quasi-linear conservative hyperbolic laws (0)-(§) 
in one-dimensional case with v = shows that for < a shock wave x = Xd{t) at the 



forward part of concentration profile Ck is formed. Simultaneously there are two fronts of 
rarefaction wave x = Xi{t) and x = Xr(t) that appear at the backward part of the profile. 
The velocities of these shock wave and the fronts of rarefaction wave are constants: 



_7i 



Vl = 71, Vr 



71 



(l + aic?)2' 



ai < 0. 



(9) 



The conductivity in the rarefaction wave is: 



lit 



1/2 



Xlit) ^ X ^ Xr(t). 



(10) 



In contrast, for the case > there is a shock wave x = Xd{t) at the backward part of the 
profile and there are two fronts of rarefaction wave x = xi{t) and x = Xr{t) on the forward 
part of the profile. The distributions shown in Fig. ^]l] do exist until the instant t = tint, 
when a direct interaction between the waves takes place. For example, for < the front 
of rarefaction wave Xrit) will catch up with the shock wave Xd{t) i.e. ^^(tint) = ^^(tint). The 
further evolution can be described analytically. Omitting details, one can notice that finally 
the concentration profile adapts a 'triangular' shape (Fig. 
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Figure 2.2: Final stage of the process 



At t ^ +00 the hight c*{t) and the base 5{t) of the 'triangle' are given (independently 
of the sign of a) as 



f M 1 

7 , 5(t)~2(«i7iMt)^/% -c*(t)5(t)~M = c?(x2-xi), (11) 



V"i7i^ 



where M is the total mass of an admixture. Notice that the 'spreading' of the concentration 
in the absence of diffusion process is: 



c*{t) = 0{t 



-l/2^ 



5{t) = 0{t 



1/2N 



t — > +CX) 



It represents a nonlinear effect of electromigration spreading that takes place due to the 
dependence of the admixture transfer velocity on the concentration. Notice that in the 
case of conventional diffusion c*{t) and 5{t) are changing similarly; that is the reason why 
experimentalists often mix up these two very different phenomena (diffusion and hyperbolic 



spreading of a concentration profile); at the same time one should take into account that 
the spreading due to diffusion is much slower than hyperbolic spreading (constants in the 
similar laws are very different!). More details on the zone evolution the zonal electrophoresis 



are given in ^ The two-dimensional version of this problem can be solved only 

numerically but we can still analytically derive that there is a shock wave of concentration 
on the forward part of the profile and a rarefaction wave of concentration at the rear part of 
a wave for a < 0. A normal coordinate to a zone boundary corresponds to the x-direction 
in one-dimensional case. 

The third factor causing the distortions of a zone reveals itself only for the domains 
with vertex points where the singularities of electric potential take place. In the simplest 
case when admixtures are absent (pure buffer solution) and a fluid is still {v = 0, = 0, 
a = (Tq), the equation divj = produces Laplace's equation for the potential Aip = with 
the boundary conditions (|^), (|^). In the vicinity of the vertex of the angle a there is a 
solution (p = Ar'^^°' cos{7r6/ a) where {r,6) are polar coordinates with the origin at a vertex. 
One can see that the gradient of such a potential possesses a singularity Vy? = 0{r'^^°'~^). In 
particular, for the domain shown in Fig. |1.1| the radial component of electric field at the 
points D and H {a = 37r/2) is very large: Ej. = 0{r~^/^). In contrary, at points A and E 
{a = vr/2) the singularities are absent and the electric field is small: Er = 0{r). Now, let all 
admixture at the initial instant be placed in the vicinity of BC and the potential difference 
(fi — (fio be such that the migration of admixture is directed towards GF. It is apparent that 
there is a fast admixture transport in the vicinities of points D and H, and a slow one in the 
vicinities of A and E. In the latter case the formation of stagnation regions is likely. This 
effect will eventually cause the strong distortion of zone shapes. 

Finally, the fourth factor is gravitational concentration convection. The difference be- 
tween the densities of an admixture and a buffer fluid generates buoyancy flows in the regions 
of inhomogeneous density. One can guess that there should be an intense fluid flow induced 
by the motion of an admixture in the vicinities of points D and H (see Fig. |1 . 1| ) . 



3 Numerical Experiments 

The problem (|l])-(0) has been solved by a direct simulation with the employing of the 
finite-difference method of markers and cells (MAC). In order to approximate the transport 
equations (^ we have used combined explicit and implicit finite-difference schemes with the 
finite differences taken in the direction opposite to a flow; the latter allows us to block mesh 
diffusion effects. To compute the flow velocity v we have used explicit schemes. Finally 
the method of 'sequence over relaxation' (SOR) with the relaxation parameter 1.37^ 1.96 
has been chosen for the computation of pressure p and potential (p (in the solving of finite- 
difference analogues of elliptic equations). It should be noticed that due to singularities in 
Vif in the vicinities of points D and H the SOR method has been essentially modified: in 
particular we introduced five subdomains with the appropriate matching conditions at their 



boundaries. In addition, the computations have also been carried out by the finite element 
method with the use of FreeFem++ and FlexPDE. The comparative relative error between 
the computations by the different methods used has been below 0.03. 

The initial concentration is prescribed in a circular zone of radius r centered at Xq, zq as 

or in a rectangular zone given by an appropriate expression. Both these distributions repre- 
sent 'almost step functions' with the parameter 6 = 100 used for the smoothing of disconti- 
nuities. 



3.1 Zone Distortion in Complex Shape Channel 

Some typical results for the channel shown in Fig. [0| are presented in Fig. p.l| ^ p.3[ In 
the case of a lighter single admixture (/i = 0.01, Pi = —0.1, 71 = —0.15, ipi — ipo = 20, 
«! = —0.4, e = 0.1) the surface levels of the concentration ci{x,z) and the streamlines 
of fiuid fiows are shown in Figs. p3.1| , p3.2| . The sequence of frames 1-6 corresponds to the 
instants t = 0.24,9.91,14.66,26.71,42.11,68.50. The potential difference chosen drives the 
zone from the line BC towards the line GF. The lengths of intervals in the channel boundary 
ABCDEFGH are: 

\BC\ =wi>0, \DE\ = 1 - wi > 0, \GF\ = ^2 > 0, \AH\ = 1 - > 0, ^1 = ^2 = ^, 

\CD\ = hi>0, \HG\ = /i2 > 0, \AB\ = 1 - > 0, \EF\ = l-hi>0, hi = h2 = -. 

3 

Combined explicit and implicit schemes have allowed us to use the mesh 60 x 60 with rather 
large steps = hy = 0.0167. It is well visible that there are strong distortions of the 
zone shape in the vicinities of points H, while in the vicinities of the points A, E the 
admixture is retarded in stagnation regions. It is interesting that the distortion is so strong 
that it causes the formation of three vortices (for the instant t = 14.66 in the frame 3, 
Fig. |3.2|) in the vicinity of the point D and the subsequent disappearing of one vortex. 




Figure 3.2: The streamfunction 
The numerical experiments have shown that the variations of the diffusion parameter e 



within the interval (0.0001 < e < 0.1) does not change the resuhs. It means that electromi- 
gration (not diffusion!) is the main reason for the zone spreading, while diffusion can affect 
the deforming of zone only in large time intervals. A series of computations have shown that 
qualitatively similar pictures with the forming of stagnation regions and the distortion of a 
zone shape have been observed in a wide interval of e. The changing of |7i||v5i — V'ol (for 
the fixed /i, jSi) leads only to the changing of a time scale for the zone passing through the 
channel. The typical size of a zone in the direction of motion (before its qualitative distor- 
tion has appeared) has been defined by the parameter ai. It agrees well with the results for 
one-dimensional case (see Figs. |2.1| , |2.2| ) at least for ai > —0.5. The zone distortion is most 
sensible to the change of the parameters hi, h2, Wi, W2 that prescribe the relative sizes of an 
electrophoretic chamber. Any quantitative description of a zone distortion is rather difficult. 
For 'symmetric' chambers with 0.05 < hi = h2 < 0.1 and 0.05 < wi = W2 < 0.1 up to 40% 
of the total mass of an admixture is trapped in the vicinity of points A and E, while for the 
case 0.3 < hi = h2 < 0.4 and 0.2 < wi = W2 < 0.3 this figure is up to 25%. The intensity of 
convection is strongly influenced by the viscosity fi or more precisely by Gr = where 
Gr is a version of Grashof 's number related to concentration. An intense (almost chaotic) 
mixing takes place for Gr > 10^. 

The frames 1 and 2, Fig. ^]3| show a stage of the separation of two admixtures that are 
heavier than the buffer (/i = 0.01, e = 0.1, ipi — ipQ = 20, Pi = 0.1, P2 = 0.1, 71 = —0.15, 
72 = —0.35, «! = —0.4, a2 = —0.4). It is interesting to see that the 'faster' admixture 
(I72I > I71I) has been undergoing the larger distortions. We have split the pictures for 
two concentrations into two frames in order to avoid visual superimposing of surface levels. 
The frame 3, Fig. p.3| shows a stage of separation of two heavier admixtures moving in 
the opposite directions (towards each other) (/i = 0.01, e = 0.1, (pi ~ ipo = 20, Pi = 0.1, 
P2 = 0.1, 7i = —0.15, 72 = 0.15, tti = —0.4, 0^2 = —0.4). It is noticeable that the motion of 
the admixture C2 against the gravity causes the greater profile distortions than the motion 
of the admixture Ci{x, z) along the gravity field. 

The described numerical experiments have shown that the motion of zones in the gravity 
field and in channels with vertex points does produce very strong zone distortions. This fact 
explains the failures of the experiments with the 'Kashtan' devices in space that were noticed 
in H, |0|. 



Figure 3.3: The surface levels of the concentrations Ci(x, z) (frame 1), C2(x, z) (frame 2), and 
Ci{x,z), C2{x,z) (frame 3) 



3.2 Zone Distortion in Rectangular Channel 

In order to demonstrate more clearly the influence of convective mixing on the zone distortion 
we have also made computations for a simple rectangular channel where the singularities of 
an electric field are absent. Constant potentials has been prescribed on side walls, while top 
and bottom walls has been insulated. The width of the channel is a = 2, while its hight is 
b = 1. The rest of parameters have been chosen as: ifi — ifQ = 20, fi = 0.01, 71 = —0.35, 
e = 0.01. Below we present the results for a single admixture with different values ai and (3i 
which show that the gravity can strongly affect the admixture transport via the buoyancy 
driven intense vortex flows near the zone. 

Figs. |3.4| , |3.5| present the surface levels for (5i = —200 (a lighter admixture) in a fluid with 
small viscosity (Gr = 2 ■ 10^) for two different values of ai at the instants t = 0.0256; 0.1024; 



0.1792; 0.2304 where one can clearly see electromigration spreading. For ai < (Fig. |3.4| ) 
the forward front exhibits the compression of surface levels (a shock wave), while the rear 
front — a rarefaction wave. For «i > (Fig. 3^) those features are opposite: the rear front 
represents a shock wave, while the forward front — a rarefaction wave. These results have 



been in agreement with the properties of motion in one-dimensional case (see Figs. 



Notice that the zone form has been essentially distorted at the initial stages of motion. The 
buoyancy effects at /5i = —200 are very weak; the admixture is slightly moving upwards 
while it is transported by an electric field. 






Figure 3.4: The surface levels, ai = —0.4 



Figure 3.5: The surface levels, ai = 0.4 



Fig. |3.6| corresponds to ai = —0.4. It shows the surface levels of concentrations at the 
instants t = 0.0512; 0.1024; 0.1280. Here we have chosen the admixture being five times 
lighter = —1000, Gr = 2 ■ 10^) than in the previous case of Figs.|3]^, ^]^; so the influence 
of the gravity here is quite essential. One can clearly see the strong interaction of the 
arising zone with the upper wall that causes the apparent splitting of a single zone into two 
(Fig. p.6|) . The next Fig. p.7| demonstrates the correspondent streamline pictures that reveal 



the forming of a vortex pair in the process of convective mixing. These results show that 
the gravity can change the zone shape drastically. 





Figure 3.6: The surface level of concentration 




Figure 3.7: The streamlines 
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